# ------------------------------------- #
# Replication code for:
#
# Rathbun, Brian, Christopher Sebastian Parker, and Caleb Pomeroy "Separate but Unequal: Ethnocentrism and Racialization 
# Explain the 'Democratic' Peace in Public Opinion," American Political Science Review.
#
# This script reproduces the Prolific survey results presented in the main text and appendices.
# ------------------------------------- #

# --- libraries --- #
library(psych)
library(ggeffects)
library(ggplot2)
library(texreg)
library(ltm)
library(plyr)
library(dplyr)
library(reshape2)
library(estimatr)
set.seed(1912)

# --- set working directory
setwd("~/Downloads/replication_files/")

# --- load data
prolific_survey <- readRDS("data/prolific_survey.rds")
nrow(prolific_survey) # N = 2659, as mentioned in the paper's text 

# --- ethnocentrism measure (from Bizumic et al)
ethno_items <- rbind(cbind(prolific_survey$eth1,
                          prolific_survey$eth2,
                          prolific_survey$eth3,
                          prolific_survey$eth4,
                          prolific_survey$eth5,
                          prolific_survey$eth6,
                          prolific_survey$eth7,
                          prolific_survey$eth8,
                          prolific_survey$eth9,
                          prolific_survey$eth10,
                          prolific_survey$eth11,
                          prolific_survey$eth12))
cronbach.alpha(ethno_items) # alpha = 0.93, as mentioned in main text
ethno_fit <- fa(ethno_items, nfactors=1, rotate="varimax", scores=TRUE, fm="pa")
ethno_fit$loadings # SS loadings 6.19, as mentioned in main text
#prolific_survey$ethno_bizumic_factor <- ethno_fit$scores[,1]

# --- prep data
prolific_survey$Y <- prolific_survey$nukes_strike
prolific_survey$W <- ifelse(prolific_survey$nukes_regime == "dem", 1, 0)
prolific_survey$Z_0 <- ifelse(prolific_survey$nukes_race == "only", 1, 0)
prolific_survey$Z_2 <- ifelse(prolific_survey$nukes_race == "nonwhite", 1, 0)
prolific_survey$E <- ifelse(prolific_survey$ethno_bizumic_factor > median(prolific_survey$ethno_bizumic_factor), 1, 0)
prolific_survey$gender <- factor(prolific_survey$gender, levels = c("nonmale", "male"))
prolific_survey$Y_threat <- prolific_survey$nukes_threat
prolific_survey$Y_immoral <- prolific_survey$nukes_immoral

# --- figure A9
{plot_demos <- prolific_survey
  plot_demos$ids <- rownames(plot_demos)
  plot_demos$age_plot <- ifelse(prolific_survey$age < 30, "18-29",
                                ifelse(prolific_survey$age >= 30 & prolific_survey$age < 40, "30-39",
                                       ifelse(prolific_survey$age >= 40 & prolific_survey$age < 50, "40-49",
                                              ifelse(prolific_survey$age >= 50 & prolific_survey$age < 60, "50-59",
                                                     ifelse(prolific_survey$age >= 60 & prolific_survey$age < 70, "60-69", "70+")))))
  
  plot_demos <- reshape2::melt(plot_demos[,c("ids", "gender", "age_plot", "ideology")], id.vars = "ids")
  plot_demos$variable <- plyr::revalue(plot_demos$variable,
                                       c("gender" = "Gender",
                                         "age_plot" = "Age",
                                         "ideology" = "Political Ideology"))
  prolific_survey$ids <- rownames(prolific_survey)
  p1 <- 
    ggplot(subset(plot_demos, variable == "Gender"), aes(x = value)) +
    geom_bar(color = "black", fill = "gray50", linewidth = .7) +
    facet_wrap(~variable, ncol = 1, scales = "free_x") +
    scale_x_discrete(labels = c("Female", "Male")) +
    theme_bw() +
    labs(x = NULL, y = "Count") +
    theme(axis.text.x = element_text(size = 13),
          axis.text.y = element_text(size = 12),
          axis.title.y = element_text(size = 13),
          strip.text.x = element_text(size = 13))
  p2 <- 
    ggplot(subset(plot_demos, variable == "Age"), aes(x = value)) +
    geom_bar(color = "black", fill = "gray50", linewidth = .7) +
    facet_wrap(~variable, ncol = 1, scales = "free_x") +
    theme_bw() +
    labs(x = NULL, y = "Count") +
    theme(axis.text.x = element_text(size = 13),
          axis.text.y = element_text(size = 12),
          axis.title.y = element_text(size = 13),
          strip.text.x = element_text(size = 13))
  p3 <- 
    ggplot(data.frame(ids = prolific_survey$ids,
                      variable = "Education",
                      value = prolific_survey$education), aes(x = value)) +
    geom_bar(color = "black", fill = "gray50", linewidth = .7) +
    facet_wrap(~variable, ncol = 1, scales = "free_x") +
    scale_x_discrete(labels = c("College\nGraduate", "Less than\nCollege")) +
    theme_bw() +
    labs(x = NULL, y = "Count") +
    theme(axis.text.x = element_text(size = 13),
          axis.text.y = element_text(size = 12),
          axis.title.y = element_text(size = 13),
          strip.text.x = element_text(size = 13))
  p4 <- 
    ggplot(subset(plot_demos, variable == "Political Ideology"), aes(x = value)) +
    geom_bar(color = "black", fill = "gray50", linewidth = .7) +
    facet_wrap(~variable, ncol = 1, scales = "free_x") +
    scale_x_discrete(labels = c("Extremely\nLiberal", "2", "3", "4", "5", "6", "Extremely\nConservative")) +
    theme_bw() +
    labs(x = NULL, y = "Count") +
    theme(axis.text.x = element_text(size = 13),
          axis.text.y = element_text(size = 12),
          axis.title.y = element_text(size = 13),
          strip.text.x = element_text(size = 13))
  
  gridExtra::grid.arrange(p1,p2,p3,p4, ncol = 1)}


# --- figure A10
{p1 <- 
  ggplot(prolific_survey, aes(ethno_bizumic_factor)) +
  geom_histogram(bins = 20, color = "black", fill ="grey90") +
  labs(x = "Ethnocentrism Factor Scores\n(Bizumic et al Superiority Scale)", y = "Count", tag = "A") +
  theme_minimal() +
  theme(legend.title = element_blank(),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 13),
        axis.title.x = element_text(size = 15, margin = margin(t=15)),
        axis.title.y = element_text(size = 14, margin = margin(r=15)),
        panel.spacing.x = unit(1.5, "lines"),
        legend.key = element_rect(color = NA, fill = NA),
        legend.spacing.y = unit(-8, "line"),
        legend.key.size = unit(1, "cm"),
        plot.tag = element_text(face = "bold", size = 22),
        plot.tag.position = c(.02, .98))
p2 <- 
  ggplot(prolific_survey, aes(ethno_bizumic_additive)) +
  geom_histogram(bins = 20, color = "black", fill ="grey90") +
  labs(x = "Ethnocentrism Additive Scores\n(Bizumic et al Superiority Scale)", y = "Count", tag = "B") +
  theme_minimal() +
  theme(legend.title = element_blank(),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 13),
        axis.title.x = element_text(size = 15, margin = margin(t=15)),
        axis.title.y = element_text(size = 14, margin = margin(r=15)),
        panel.spacing.x = unit(1.5, "lines"),
        legend.key = element_rect(color = NA, fill = NA),
        legend.spacing.y = unit(-8, "line"),
        legend.key.size = unit(1, "cm"),
        plot.tag = element_text(face = "bold", size = 22),
        plot.tag.position = c(.02, .98))

gridExtra::grid.arrange(p1, p2, nrow = 1)}


# --- prolific analysis: full sample ---#
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95 <- lm_robust(Y ~ W, data = prolific_survey, 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95 <- lm_robust(Y ~ W * Z_0, data = prolific_survey, 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
#1-(-0.30576 / -0.5568) # effect is 45.1% smaller when country is nonwhite, as mentioned in the paper's text

# --- table A11
screenreg(list(lm_dem_95, lm_nonwhite_95))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90 <- lm_robust(Y ~ W, data = prolific_survey, 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90 <- lm_robust(Y ~ W * Z_0, data = prolific_survey, 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from the full sample at alpha .05 and .10 levels; store in DF to plot later.
plot_df_full_sample <- 
  data.frame(effect_type = c("ATE", "ACDE", "Eliminated Effect"),
             est = c(lm_dem_95$coefficients[2], lm_nonwhite_95$coefficients[2], lm_nonwhite_95$coefficients[4]),
             ci_95_high = c(lm_dem_95$conf.high[2], lm_nonwhite_95$conf.high[2], lm_nonwhite_95$conf.high[4]),
             ci_95_low = c(lm_dem_95$conf.low[2], lm_nonwhite_95$conf.low[2], lm_nonwhite_95$conf.low[4]),
             ci_90_high = c(lm_dem_90$conf.high[2], lm_nonwhite_90$conf.high[2], lm_nonwhite_90$conf.high[4]),
             ci_90_low = c(lm_dem_90$conf.low[2], lm_nonwhite_90$conf.low[2], lm_nonwhite_90$conf.low[4]),
             country_race = c("Country Race Unspecified","Country Race = Nonwhite",  "Country Race = Nonwhite"),
             type = "full_sample")


# --- prolific analysis: high ethnocentrism results ---#
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_he <- lm_robust(Y ~ W, data = subset(prolific_survey, E == 1), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_he <- lm_robust(Y ~ W * Z_0, data = subset(prolific_survey, E == 1), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table A12
screenreg(list(lm_dem_95_he, lm_nonwhite_95_he))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_he <- lm_robust(Y ~ W, data = subset(prolific_survey, E == 1), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_he <- lm_robust(Y ~ W * Z_0, data = subset(prolific_survey, E == 1), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from the high ethno respondents at alpha .05 and .10 levels; store in DF to plot later.
plot_df_high_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE","Eliminated Effect"),
             est = c(lm_dem_95_he$coefficients[2], lm_nonwhite_95_he$coefficients[2], lm_nonwhite_95_he$coefficients[4]),
             ci_95_high = c(lm_dem_95_he$conf.high[2], lm_nonwhite_95_he$conf.high[2],  lm_nonwhite_95_he$conf.high[4]),
             ci_95_low = c(lm_dem_95_he$conf.low[2], lm_nonwhite_95_he$conf.low[2],  lm_nonwhite_95_he$conf.low[4]),
             ci_90_high = c(lm_dem_90_he$conf.high[2],  lm_nonwhite_90_he$conf.high[2], lm_nonwhite_90_he$conf.high[4]),
             ci_90_low = c(lm_dem_90_he$conf.low[2],  lm_nonwhite_90_he$conf.low[2],  lm_nonwhite_90_he$conf.low[4]),
             country_race = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = Nonwhite"),
             type = "high_ethno")

# --- vary the ethnocentrism threshold
# in the main text, we present results with ethnocentrism binned at the median
# here, we simply note that this result is robust beyond this median split
# nonwhite information significantly eliminates the democratic peace effect for approximately 76% of the sample, as mentioned in the paper's text

screenreg(lm_robust(Y ~ W * Z_0, data = subset(prolific_survey, ethno_bizumic_factor >= quantile(prolific_survey$ethno_bizumic_factor, c(.24))), 
                    subset = (Z_0 == 1 | Z_2 == 1),
                    alpha = .05))

# --- prolific analysis: low ethnocentrism results ---#
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_le <- lm_robust(Y ~ W, data = subset(prolific_survey, E == 0), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_le <- lm_robust(Y ~ W * Z_0, data = subset(prolific_survey, E == 0), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))

#1-(-0.33629 / -0.4013) # among respondents lower in ethnocentrism, effect is 16.2% smaller when country is nonwhite, as mentioned in the paper's text

# --- table A13
screenreg(list(lm_dem_95_le, lm_nonwhite_95_le))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_le <- lm_robust(Y ~ W, data = subset(prolific_survey, E == 0), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_le <- lm_robust(Y ~ W * Z_0, data = subset(prolific_survey, E == 0), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from the low ethno respondents at alpha .05 and .10 levels; store in DF to plot later.
plot_df_low_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE","Eliminated Effect"),
             est = c(lm_dem_95_le$coefficients[2], lm_nonwhite_95_le$coefficients[2], lm_nonwhite_95_le$coefficients[4]),
             ci_95_high = c(lm_dem_95_le$conf.high[2], lm_nonwhite_95_le$conf.high[2],  lm_nonwhite_95_le$conf.high[4]),
             ci_95_low = c(lm_dem_95_le$conf.low[2], lm_nonwhite_95_le$conf.low[2],  lm_nonwhite_95_le$conf.low[4]),
             ci_90_high = c(lm_dem_90_le$conf.high[2],  lm_nonwhite_90_le$conf.high[2], lm_nonwhite_90_le$conf.high[4]),
             ci_90_low = c(lm_dem_90_le$conf.low[2],  lm_nonwhite_90_le$conf.low[2],  lm_nonwhite_90_le$conf.low[4]),
             country_race = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = Nonwhite"),
             type = "low_ethno")

# combine the three sets of results
{plot_df <- rbind(plot_df_full_sample,
                 plot_df_high_ethno,
                 plot_df_low_ethno)

plot_df$effect_type <- factor(plot_df$effect_type, levels = c("Eliminated Effect", "ACDE", "ATE"))
levels(plot_df$effect_type) <- c("Effect of Democracy\nEliminated by Race", "Effect of Democracy\nFixing Race", "Effect of Democracy")
plot_df$country_race <- factor(plot_df$country_race, levels = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = White"))
plot_df$type <- factor(plot_df$type, levels = c("full_sample","high_ethno", "low_ethno"))
levels(plot_df$type) <- c("All Respondents\n(Prolific Sample)", "Respondents Above\nEthnocentrism Median",
                          "Respondents Below\nEthnocentrism Median")

# --- figure 3 from the main text
ggplot(plot_df, aes(x = est, y = effect_type, group = country_race, color = country_race)) +
  geom_vline(xintercept = 0, size = .8, color = "black") +
  geom_linerange(aes(xmin = ci_95_low, xmax = ci_95_high), size = 1.3, color = "gray10", position = position_dodge(width = .5)) +
  geom_linerange(aes(xmin = ci_90_low, xmax = ci_90_high), size = 2.3, position = position_dodge(width = .5)) +
  geom_point(size = 5.5, color = "black", position = position_dodge(width = .5)) +
  geom_point(size = 4.5,position = position_dodge(width = .5)) +  
  facet_wrap(~type) +
  scale_color_manual(values=c("gray30", "#92978a", "#c5bec2")) +
  labs(y = NULL, x = "Estimate") +
  theme_minimal() +
  theme(legend.title = element_blank(),
        legend.text = element_text(size = 16),
        axis.text.y = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.line.y.right = element_blank(),
        axis.title.x = element_text(size = 22, margin = margin(t=15)),
        panel.spacing.x = unit(1.5, "lines"),
        panel.grid.minor.x = element_blank(),
        strip.text.x = element_text(size = 18, margin = margin (13, 0, 15, 0)),
        strip.background = element_rect(color = "black", size = 1.5, fill = "#f3f3f1"),
        legend.key = element_rect(color = NA, fill = NA),
        legend.spacing.y = unit(-8, "line"),
        legend.key.size = unit(1, "cm"))}

# --- prolific analysis: full sample (with covariates) ---#
# here, we simply re-run the above models with the addition of right hand side covariates
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95 <- lm_robust(Y ~ W + gender + race + pid + ideology + education + age, data = prolific_survey, 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95 <- lm_robust(Y ~ W * Z_0 + gender + race + pid + ideology + education + age, data = prolific_survey, 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))

# --- table A14
screenreg(list(lm_dem_95, lm_nonwhite_95))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90 <- lm_robust(Y ~ W + gender + race + pid + ideology + education + age, data = prolific_survey, 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90 <- lm_robust(Y ~ W * Z_0 + gender + race + pid + ideology + education + age, data = prolific_survey, 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# again, results from the full sample at alpha .05 and .10 levels; this time w/ covariates; store in DF to plot later.
plot_df_full_sample <- 
  data.frame(effect_type = c("ATE", "ACDE", "Eliminated Effect"),
             est = c(lm_dem_95$coefficients[2], lm_nonwhite_95$coefficients[2], lm_nonwhite_95$coefficients[length(lm_nonwhite_95$coefficients)]),
             ci_95_high = c(lm_dem_95$conf.high[2], lm_nonwhite_95$conf.high[2], lm_nonwhite_95$conf.high[length(lm_nonwhite_95$coefficients)]),
             ci_95_low = c(lm_dem_95$conf.low[2], lm_nonwhite_95$conf.low[2], lm_nonwhite_95$conf.low[length(lm_nonwhite_95$coefficients)]),
             ci_90_high = c(lm_dem_90$conf.high[2], lm_nonwhite_90$conf.high[2], lm_nonwhite_90$conf.high[length(lm_nonwhite_90$coefficients)]),
             ci_90_low = c(lm_dem_90$conf.low[2], lm_nonwhite_90$conf.low[2], lm_nonwhite_90$conf.low[length(lm_nonwhite_90$coefficients)]),
             country_race = c("Country Race Unspecified","Country Race = Nonwhite",  "Country Race = Nonwhite"),
             type = "full_sample")

# --- prolific analysis: high ethnocentrism results ---#
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_he <- lm_robust(Y ~ W + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 1), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_he <- lm_robust(Y ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 1), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table A15
screenreg(list(lm_dem_95_he, lm_nonwhite_95_he))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_he <- lm_robust(Y ~ W + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 1), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_he <- lm_robust(Y ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 1), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# again, results from the high ethno respondents at alpha .05 and .10 levels; this time w/ covariates; store in DF to plot later.
plot_df_high_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE","Eliminated Effect"),
             est = c(lm_dem_95_he$coefficients[2], lm_nonwhite_95_he$coefficients[2], lm_nonwhite_95_he$coefficients[length(lm_nonwhite_95_he$coefficients)]),
             ci_95_high = c(lm_dem_95_he$conf.high[2], lm_nonwhite_95_he$conf.high[2],  lm_nonwhite_95_he$conf.high[length(lm_nonwhite_95_he$coefficients)]),
             ci_95_low = c(lm_dem_95_he$conf.low[2], lm_nonwhite_95_he$conf.low[2],  lm_nonwhite_95_he$conf.low[length(lm_nonwhite_95_he$coefficients)]),
             ci_90_high = c(lm_dem_90_he$conf.high[2],  lm_nonwhite_90_he$conf.high[2], lm_nonwhite_90_he$conf.high[length(lm_nonwhite_90_he$coefficients)]),
             ci_90_low = c(lm_dem_90_he$conf.low[2],  lm_nonwhite_90_he$conf.low[2],  lm_nonwhite_90_he$conf.low[length(lm_nonwhite_90_he$coefficients)]),
             country_race = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = Nonwhite"),
             type = "high_ethno")

# --- prolific analysis: low ethnocentrism results ---#
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_le <- lm_robust(Y ~ W + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 0), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_le <- lm_robust(Y ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 0), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table A16
screenreg(list(lm_dem_95_le, lm_nonwhite_95_le))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_le <- lm_robust(Y ~ W + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 0), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_le <- lm_robust(Y ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 0), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# again, results from the low ethno respondents at alpha .05 and .10 levels; this time w/ covariates; store in DF to plot later.
plot_df_low_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE","Eliminated Effect"),
             est = c(lm_dem_95_le$coefficients[2], lm_nonwhite_95_le$coefficients[2], lm_nonwhite_95_le$coefficients[length(lm_nonwhite_95_le$coefficients)]),
             ci_95_high = c(lm_dem_95_le$conf.high[2], lm_nonwhite_95_le$conf.high[2],  lm_nonwhite_95_le$conf.high[length(lm_nonwhite_95_le$coefficients)]),
             ci_95_low = c(lm_dem_95_le$conf.low[2], lm_nonwhite_95_le$conf.low[2],  lm_nonwhite_95_le$conf.low[length(lm_nonwhite_95_le$coefficients)]),
             ci_90_high = c(lm_dem_90_le$conf.high[2],  lm_nonwhite_90_le$conf.high[2], lm_nonwhite_90_le$conf.high[length(lm_nonwhite_90_le$coefficients)]),
             ci_90_low = c(lm_dem_90_le$conf.low[2],  lm_nonwhite_90_le$conf.low[2],  lm_nonwhite_90_le$conf.low[length(lm_nonwhite_90_le$coefficients)]),
             country_race = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = Nonwhite"),
             type = "low_ethno")

# combine the three sets of results w/ covariates modeled
{plot_df <- rbind(plot_df_full_sample,
                 plot_df_high_ethno,
                 plot_df_low_ethno)

plot_df$effect_type <- factor(plot_df$effect_type, levels = c("Eliminated Effect", "ACDE", "ATE"))
levels(plot_df$effect_type) <- c("Effect of Democracy\nEliminated by Race", "Effect of Democracy\nFixing Race", "Effect of Democracy")
plot_df$country_race <- factor(plot_df$country_race, levels = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = White"))
plot_df$type <- factor(plot_df$type, levels = c("full_sample","high_ethno", "low_ethno"))
levels(plot_df$type) <- c("All Respondents\n(Prolific Sample)", "Respondents Above\nEthnocentrism Median",
                          "Respondents Below\nEthnocentrism Median")
# --- figure A11
ggplot(plot_df, aes(x = est, y = effect_type, group = country_race, color = country_race)) +
  geom_vline(xintercept = 0, size = .8, color = "black") +
  geom_linerange(aes(xmin = ci_95_low, xmax = ci_95_high), size = 1.3, color = "gray10", position = position_dodge(width = .5)) +
  geom_linerange(aes(xmin = ci_90_low, xmax = ci_90_high), size = 2.3, position = position_dodge(width = .5)) +
  geom_point(size = 5.5, color = "black", position = position_dodge(width = .5)) +
  geom_point(size = 4.5,position = position_dodge(width = .5)) +  
  facet_wrap(~type) +
  scale_color_manual(values=c("gray30", "#92978a", "#c5bec2")) +
  labs(y = NULL, x = "Estimate") +
  theme_minimal() +
  theme(legend.title = element_blank(),
        legend.text = element_text(size = 16),
        axis.text.y = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.line.y.right = element_blank(),
        axis.title.x = element_text(size = 22, margin = margin(t=15)),
        panel.spacing.x = unit(1.5, "lines"),
        panel.grid.minor.x = element_blank(),
        strip.text.x = element_text(size = 18, margin = margin (13, 0, 15, 0)),
        strip.background = element_rect(color = "black", size = 1.5, fill = "#f3f3f1"),
        legend.key = element_rect(color = NA, fill = NA),
        legend.spacing.y = unit(-8, "line"),
        legend.key.size = unit(1, "cm"))}

# --- prolific analysis: threat perception and immorality mechanisms ---#
# same estimates as above, this time for the threat and immorality mechanisms identified in Tomz and Weeks;
# --- starting with threat perception estimates, full sample
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95 <- lm_robust(Y_threat ~ W + gender + race + pid + ideology + education + age, data = prolific_survey, 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95 <- lm_robust(Y_threat ~ W * Z_0 + gender + race + pid + ideology + education + age, data = prolific_survey, 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table B7
screenreg(list(lm_dem_95, lm_nonwhite_95))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90 <- lm_robust(Y_threat ~ W + gender + race + pid + ideology + education + age, data = prolific_survey, 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90 <- lm_robust(Y_threat ~ W * Z_0 + gender + race + pid + ideology + education + age, data = prolific_survey, 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))
# results from the full sample at alpha .05 and .10 levels; this time for threat perception DV; store in DF to plot later.
plot_df_full_sample <- 
  data.frame(effect_type = c("ATE", "ACDE", "Eliminated Effect"),
             est = c(lm_dem_95$coefficients[2], lm_nonwhite_95$coefficients[2], lm_nonwhite_95$coefficients[length(lm_nonwhite_95$coefficients)]),
             ci_95_high = c(lm_dem_95$conf.high[2], lm_nonwhite_95$conf.high[2], lm_nonwhite_95$conf.high[length(lm_nonwhite_95$coefficients)]),
             ci_95_low = c(lm_dem_95$conf.low[2], lm_nonwhite_95$conf.low[2], lm_nonwhite_95$conf.low[length(lm_nonwhite_95$coefficients)]),
             ci_90_high = c(lm_dem_90$conf.high[2], lm_nonwhite_90$conf.high[2], lm_nonwhite_90$conf.high[length(lm_nonwhite_90$coefficients)]),
             ci_90_low = c(lm_dem_90$conf.low[2], lm_nonwhite_90$conf.low[2], lm_nonwhite_90$conf.low[length(lm_nonwhite_90$coefficients)]),
             country_race = c("Country Race Unspecified","Country Race = Nonwhite",  "Country Race = Nonwhite"),
             type = "full_sample")


# --- threat perception estimates, high ethno respondents
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_he <- lm_robust(Y_threat ~ W + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 1), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_he <- lm_robust(Y_threat ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 1), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table B8
screenreg(list(lm_dem_95_he, lm_nonwhite_95_he))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_he <- lm_robust(Y_threat ~ W + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 1), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_he <- lm_robust(Y_threat ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 1), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from the high ethno respondents at alpha .05 and .10 levels; this time for threat perception DV; store in DF to plot later.
plot_df_high_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE","Eliminated Effect"),
             est = c(lm_dem_95_he$coefficients[2], lm_nonwhite_95_he$coefficients[2], lm_nonwhite_95_he$coefficients[length(lm_nonwhite_95_he$coefficients)]),
             ci_95_high = c(lm_dem_95_he$conf.high[2], lm_nonwhite_95_he$conf.high[2],  lm_nonwhite_95_he$conf.high[length(lm_nonwhite_95_he$coefficients)]),
             ci_95_low = c(lm_dem_95_he$conf.low[2], lm_nonwhite_95_he$conf.low[2],  lm_nonwhite_95_he$conf.low[length(lm_nonwhite_95_he$coefficients)]),
             ci_90_high = c(lm_dem_90_he$conf.high[2],  lm_nonwhite_90_he$conf.high[2], lm_nonwhite_90_he$conf.high[length(lm_nonwhite_90_he$coefficients)]),
             ci_90_low = c(lm_dem_90_he$conf.low[2],  lm_nonwhite_90_he$conf.low[2],  lm_nonwhite_90_he$conf.low[length(lm_nonwhite_90_he$coefficients)]),
             country_race = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = Nonwhite"),
             type = "high_ethno")

# --- threat perception estimates, low ethno respondents
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_le <- lm_robust(Y_threat ~ W + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 0), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_le <- lm_robust(Y_threat ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 0), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))

# --- table B9
screenreg(list(lm_dem_95_le, lm_nonwhite_95_le))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_le <- lm_robust(Y_threat ~ W + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 0), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_le <- lm_robust(Y_threat ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 0), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from the low ethno respondents at alpha .05 and .10 levels; this time for threat perception DV; store in DF to plot later.
plot_df_low_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE","Eliminated Effect"),
             est = c(lm_dem_95_le$coefficients[2], lm_nonwhite_95_le$coefficients[2], lm_nonwhite_95_le$coefficients[length(lm_nonwhite_95_le$coefficients)]),
             ci_95_high = c(lm_dem_95_le$conf.high[2], lm_nonwhite_95_le$conf.high[2],  lm_nonwhite_95_le$conf.high[length(lm_nonwhite_95_le$coefficients)]),
             ci_95_low = c(lm_dem_95_le$conf.low[2], lm_nonwhite_95_le$conf.low[2],  lm_nonwhite_95_le$conf.low[length(lm_nonwhite_95_le$coefficients)]),
             ci_90_high = c(lm_dem_90_le$conf.high[2],  lm_nonwhite_90_le$conf.high[2], lm_nonwhite_90_le$conf.high[length(lm_nonwhite_90_le$coefficients)]),
             ci_90_low = c(lm_dem_90_le$conf.low[2],  lm_nonwhite_90_le$conf.low[2],  lm_nonwhite_90_le$conf.low[length(lm_nonwhite_90_le$coefficients)]),
             country_race = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = Nonwhite"),
             type = "low_ethno")

# combine the three sets of threat perception results
{plot_df <- rbind(plot_df_full_sample,
                 plot_df_high_ethno,
                 plot_df_low_ethno)

plot_df$effect_type <- factor(plot_df$effect_type, levels = c("Eliminated Effect", "ACDE", "ATE"))
levels(plot_df$effect_type) <- c("Effect of Democracy\nEliminated by Race", "Effect of Democracy\nFixing Race", "Effect of Democracy")
plot_df$country_race <- factor(plot_df$country_race, levels = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = White"))
plot_df$type <- factor(plot_df$type, levels = c("full_sample","high_ethno", "low_ethno"))
levels(plot_df$type) <- c("All Respondents\n(Prolific Sample)", "Respondents Above\nEthnocentrism Median",
                          "Respondents Below\nEthnocentrism Median")
# --- figure A12
ggplot(plot_df, aes(x = est, y = effect_type, group = country_race, color = country_race)) +
  geom_vline(xintercept = 0, size = .8, color = "black") +
  geom_linerange(aes(xmin = ci_95_low, xmax = ci_95_high), size = 1.3, color = "gray10", position = position_dodge(width = .5)) +
  geom_linerange(aes(xmin = ci_90_low, xmax = ci_90_high), size = 2.3, position = position_dodge(width = .5)) +
  geom_point(size = 5.5, color = "black", position = position_dodge(width = .5)) +
  geom_point(size = 4.5,position = position_dodge(width = .5)) +  
  facet_wrap(~type) +
  scale_color_manual(values=c("gray30", "#92978a", "#c5bec2")) +
  labs(y = NULL, x = "Threat Perception Estimate") +
  theme_minimal() +
  theme(legend.title = element_blank(),
        legend.text = element_text(size = 16),
        axis.text.y = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.line.y.right = element_blank(),
        axis.title.x = element_text(size = 22, margin = margin(t=15)),
        panel.spacing.x = unit(1.5, "lines"),
        panel.grid.minor.x = element_blank(),
        strip.text.x = element_text(size = 18, margin = margin (13, 0, 15, 0)),
        strip.background = element_rect(color = "black", size = 1.5, fill = "#f3f3f1"),
        legend.key = element_rect(color = NA, fill = NA),
        legend.spacing.y = unit(-8, "line"),
        legend.key.size = unit(1, "cm"))}


# --- next, immorality estimates, full sample
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95 <- lm_robust(Y_immoral ~ W + gender + race + pid + ideology + education + age, data = prolific_survey, 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95 <- lm_robust(Y_immoral ~ W * Z_0 + gender + race + pid + ideology + education + age, data = prolific_survey, 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))

# --- table B10
screenreg(list(lm_dem_95, lm_nonwhite_95))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90 <- lm_robust(Y_immoral ~ W + gender + race + pid + ideology + education + age, data = prolific_survey, 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90 <- lm_robust(Y_immoral ~ W * Z_0 + gender + race + pid + ideology + education + age, data = prolific_survey, 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from the full sample at alpha .05 and .10 levels; this time for immorality DV; store in DF to plot later.
plot_df_full_sample <- 
  data.frame(effect_type = c("ATE", "ACDE", "Eliminated Effect"),
             est = c(lm_dem_95$coefficients[2], lm_nonwhite_95$coefficients[2], lm_nonwhite_95$coefficients[length(lm_nonwhite_95$coefficients)]),
             ci_95_high = c(lm_dem_95$conf.high[2], lm_nonwhite_95$conf.high[2], lm_nonwhite_95$conf.high[length(lm_nonwhite_95$coefficients)]),
             ci_95_low = c(lm_dem_95$conf.low[2], lm_nonwhite_95$conf.low[2], lm_nonwhite_95$conf.low[length(lm_nonwhite_95$coefficients)]),
             ci_90_high = c(lm_dem_90$conf.high[2], lm_nonwhite_90$conf.high[2], lm_nonwhite_90$conf.high[length(lm_nonwhite_90$coefficients)]),
             ci_90_low = c(lm_dem_90$conf.low[2], lm_nonwhite_90$conf.low[2], lm_nonwhite_90$conf.low[length(lm_nonwhite_90$coefficients)]),
             country_race = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = Nonwhite"),
             type = "full_sample")


# --- immorality estimates, high ethno respondents
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_he <- lm_robust(Y_immoral ~ W + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 1), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_he <- lm_robust(Y_immoral ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 1), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table B11
screenreg(list(lm_dem_95_he, lm_nonwhite_95_he))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_he <- lm_robust(Y_immoral ~ W + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 1), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_he <- lm_robust(Y_immoral ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 1), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results for high ethno respondents at alpha .05 and .10 levels; this time for immorality DV; store in DF to plot later.
plot_df_high_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE","Eliminated Effect"),
             est = c(lm_dem_95_he$coefficients[2], lm_nonwhite_95_he$coefficients[2], lm_nonwhite_95_he$coefficients[length(lm_nonwhite_95_he$coefficients)]),
             ci_95_high = c(lm_dem_95_he$conf.high[2], lm_nonwhite_95_he$conf.high[2],  lm_nonwhite_95_he$conf.high[length(lm_nonwhite_95_he$coefficients)]),
             ci_95_low = c(lm_dem_95_he$conf.low[2], lm_nonwhite_95_he$conf.low[2],  lm_nonwhite_95_he$conf.low[length(lm_nonwhite_95_he$coefficients)]),
             ci_90_high = c(lm_dem_90_he$conf.high[2],  lm_nonwhite_90_he$conf.high[2], lm_nonwhite_90_he$conf.high[length(lm_nonwhite_90_he$coefficients)]),
             ci_90_low = c(lm_dem_90_he$conf.low[2],  lm_nonwhite_90_he$conf.low[2],  lm_nonwhite_90_he$conf.low[length(lm_nonwhite_90_he$coefficients)]),
             country_race = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = Nonwhite"),
             type = "high_ethno")

# --- immorality estimates, low ethno respondents
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_le <- lm_robust(Y_immoral ~ W + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 0), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_le <- lm_robust(Y_immoral ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 0), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))

# --- table B12
screenreg(list(lm_dem_95_le, lm_nonwhite_95_le))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_le <- lm_robust(Y_immoral ~ W + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 0), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_le <- lm_robust(Y_immoral ~ W * Z_0 + gender + race + pid + ideology + education + age, data = subset(prolific_survey, E == 0), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results for low ethno respondents at alpha .05 and .10 levels; this time for immorality DV; store in DF to plot later.
plot_df_low_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE","Eliminated Effect"),
             est = c(lm_dem_95_le$coefficients[2], lm_nonwhite_95_le$coefficients[2], lm_nonwhite_95_le$coefficients[length(lm_nonwhite_95_le$coefficients)]),
             ci_95_high = c(lm_dem_95_le$conf.high[2], lm_nonwhite_95_le$conf.high[2],  lm_nonwhite_95_le$conf.high[length(lm_nonwhite_95_le$coefficients)]),
             ci_95_low = c(lm_dem_95_le$conf.low[2], lm_nonwhite_95_le$conf.low[2],  lm_nonwhite_95_le$conf.low[length(lm_nonwhite_95_le$coefficients)]),
             ci_90_high = c(lm_dem_90_le$conf.high[2],  lm_nonwhite_90_le$conf.high[2], lm_nonwhite_90_le$conf.high[length(lm_nonwhite_90_le$coefficients)]),
             ci_90_low = c(lm_dem_90_le$conf.low[2],  lm_nonwhite_90_le$conf.low[2],  lm_nonwhite_90_le$conf.low[length(lm_nonwhite_90_le$coefficients)]),
             country_race = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = Nonwhite"),
             type = "low_ethno")

# combine the three sets of immorality results
{plot_df <- rbind(plot_df_full_sample,
                 plot_df_high_ethno,
                 plot_df_low_ethno)

plot_df$effect_type <- factor(plot_df$effect_type, levels = c("Eliminated Effect", "ACDE", "ATE"))
levels(plot_df$effect_type) <- c("Effect of Democracy\nEliminated by Race", "Effect of Democracy\nFixing Race", "Effect of Democracy")
plot_df$country_race <- factor(plot_df$country_race, levels = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = White"))
plot_df$type <- factor(plot_df$type, levels = c("full_sample","high_ethno", "low_ethno"))
levels(plot_df$type) <- c("All Respondents\n(Prolific Sample)", "Respondents Above\nEthnocentrism Median",
                          "Respondents Below\nEthnocentrism Median")
# --- figure A13
ggplot(plot_df, aes(x = est, y = effect_type, group = country_race, color = country_race)) +
  geom_vline(xintercept = 0, size = .8, color = "black") +
  geom_linerange(aes(xmin = ci_95_low, xmax = ci_95_high), size = 1.3, color = "gray10", position = position_dodge(width = .5)) +
  geom_linerange(aes(xmin = ci_90_low, xmax = ci_90_high), size = 2.3, position = position_dodge(width = .5)) +
  geom_point(size = 5.5, color = "black", position = position_dodge(width = .5)) +
  geom_point(size = 4.5,position = position_dodge(width = .5)) +  
  facet_wrap(~type) +
  scale_color_manual(values=c("gray30", "#92978a", "#c5bec2")) +
  labs(y = NULL, x = "Perceived Immorality Estimate") +
  theme_minimal() +
  theme(legend.title = element_blank(),
        legend.text = element_text(size = 16),
        axis.text.y = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.line.y.right = element_blank(),
        axis.title.x = element_text(size = 22, margin = margin(t=15)),
        panel.spacing.x = unit(1.5, "lines"),
        panel.grid.minor.x = element_blank(),
        strip.text.x = element_text(size = 18, margin = margin (13, 0, 15, 0)),
        strip.background = element_rect(color = "black", size = 1.5, fill = "#f3f3f1"),
        legend.key = element_rect(color = NA, fill = NA),
        legend.spacing.y = unit(-8, "line"),
        legend.key.size = unit(1, "cm"))}


# --- prolific analysis: main effect of race ---#
# as noted in the dataverse appendix, we find little evidence of a main effect of race on support for a strike
t.test(subset(prolific_survey, nukes_race == "nonwhite")$Y,
       subset(prolific_survey, nukes_race == "only")$Y)

# --- prolific analysis: nonwhite identifying respondents ---#
# here, we re-estimate the primary experimental results for the nonwhite identifying respondents in our sample
# however, note: the experiment was not designed with sufficient statistical power to identify these effects.
# --- all nonwhite identifying respondents
nrow(subset(prolific_survey, race == "nonwhite"))
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95 <- lm_robust(Y ~ W, data = subset(prolific_survey, race == "nonwhite"), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95 <- lm_robust(Y ~ W * Z_0, data = subset(prolific_survey, race == "nonwhite"), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table B4
screenreg(list(lm_dem_95, lm_nonwhite_95))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90 <- lm_robust(Y ~ W, data = subset(prolific_survey, race == "nonwhite"), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90 <- lm_robust(Y ~ W * Z_0, data = subset(prolific_survey, race == "nonwhite"), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from nonwhite identifying respondents at alpha .05 and .10 levels; store in DF to plot later.
plot_df_full_sample <- 
  data.frame(effect_type = c("ATE", "ACDE", "Eliminated Effect"),
             est = c(lm_dem_95$coefficients[2], lm_nonwhite_95$coefficients[2], lm_nonwhite_95$coefficients[4]),
             ci_95_high = c(lm_dem_95$conf.high[2], lm_nonwhite_95$conf.high[2], lm_nonwhite_95$conf.high[4]),
             ci_95_low = c(lm_dem_95$conf.low[2], lm_nonwhite_95$conf.low[2], lm_nonwhite_95$conf.low[4]),
             ci_90_high = c(lm_dem_90$conf.high[2], lm_nonwhite_90$conf.high[2], lm_nonwhite_90$conf.high[4]),
             ci_90_low = c(lm_dem_90$conf.low[2], lm_nonwhite_90$conf.low[2], lm_nonwhite_90$conf.low[4]),
             country_race = c("Country Race Unspecified","Country Race = Nonwhite",  "Country Race = Nonwhite"),
             type = "full_sample")

# --- nonwhite identifying respondents high in ethnocentrism
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_he <- lm_robust(Y ~ W, data = subset(prolific_survey, E == 1 & race == "nonwhite"), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_he <- lm_robust(Y ~ W * Z_0, data = subset(prolific_survey, E == 1 & race == "nonwhite"), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table B5
screenreg(list(lm_dem_95_he, lm_nonwhite_95_he))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_he <- lm_robust(Y ~ W, data = subset(prolific_survey, E == 1 & race == "nonwhite"), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_he <- lm_robust(Y ~ W * Z_0, data = subset(prolific_survey, E == 1 & race == "nonwhite"), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from nonwhite identifying respondents above ethno median, at alpha .05 and .10 levels; store in DF to plot later.
plot_df_high_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE","Eliminated Effect"),
             est = c(lm_dem_95_he$coefficients[2], lm_nonwhite_95_he$coefficients[2], lm_nonwhite_95_he$coefficients[4]),
             ci_95_high = c(lm_dem_95_he$conf.high[2], lm_nonwhite_95_he$conf.high[2],  lm_nonwhite_95_he$conf.high[4]),
             ci_95_low = c(lm_dem_95_he$conf.low[2], lm_nonwhite_95_he$conf.low[2],  lm_nonwhite_95_he$conf.low[4]),
             ci_90_high = c(lm_dem_90_he$conf.high[2],  lm_nonwhite_90_he$conf.high[2], lm_nonwhite_90_he$conf.high[4]),
             ci_90_low = c(lm_dem_90_he$conf.low[2],  lm_nonwhite_90_he$conf.low[2],  lm_nonwhite_90_he$conf.low[4]),
             country_race = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = Nonwhite"),
             type = "high_ethno")

# --- nonwhite identifying respondents low in ethnocentrism
# --- natural mediator arm (alpha = .05)
summary(lm_dem_95_le <- lm_robust(Y ~ W, data = subset(prolific_survey, E == 0 & race == "nonwhite"), 
                               subset = Z_0 == 1,
                               alpha = .05))
# --- nonwhite mediator arm (alpha = .05)
summary(lm_nonwhite_95_le <- lm_robust(Y ~ W * Z_0, data = subset(prolific_survey, E == 0 & race == "nonwhite"), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .05))
# --- table B6
screenreg(list(lm_dem_95_le, lm_nonwhite_95_le))

# --- natural mediator arm (alpha = .10)
summary(lm_dem_90_le <- lm_robust(Y ~ W, data = subset(prolific_survey, E == 0 & race == "nonwhite"), 
                               subset = Z_0 == 1,
                               alpha = .10))
# --- nonwhite mediator arm (alpha = .10)
summary(lm_nonwhite_90_le <- lm_robust(Y ~ W * Z_0, data = subset(prolific_survey, E == 0 & race == "nonwhite"), 
                                    subset = (Z_0 == 1 | Z_2 == 1),
                                    alpha = .10))

# results from nonwhite identifying respondents below ethno median, at alpha .05 and .10 levels; store in DF to plot later.
plot_df_low_ethno <- 
  data.frame(effect_type = c("ATE", "ACDE","Eliminated Effect"),
             est = c(lm_dem_95_le$coefficients[2], lm_nonwhite_95_le$coefficients[2], lm_nonwhite_95_le$coefficients[4]),
             ci_95_high = c(lm_dem_95_le$conf.high[2], lm_nonwhite_95_le$conf.high[2],  lm_nonwhite_95_le$conf.high[4]),
             ci_95_low = c(lm_dem_95_le$conf.low[2], lm_nonwhite_95_le$conf.low[2],  lm_nonwhite_95_le$conf.low[4]),
             ci_90_high = c(lm_dem_90_le$conf.high[2],  lm_nonwhite_90_le$conf.high[2], lm_nonwhite_90_le$conf.high[4]),
             ci_90_low = c(lm_dem_90_le$conf.low[2],  lm_nonwhite_90_le$conf.low[2],  lm_nonwhite_90_le$conf.low[4]),
             country_race = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = Nonwhite"),
             type = "low_ethno")

# combine the three sets of results
{plot_df <- rbind(plot_df_full_sample,
                 plot_df_high_ethno,
                 plot_df_low_ethno)

plot_df$effect_type <- factor(plot_df$effect_type, levels = c("Eliminated Effect", "ACDE", "ATE"))
levels(plot_df$effect_type) <- c("Effect of Democracy\nEliminated by Race", "Effect of Democracy\nFixing Race", "Effect of Democracy")
plot_df$country_race <- factor(plot_df$country_race, levels = c("Country Race Unspecified", "Country Race = Nonwhite", "Country Race = White"))
plot_df$type <- factor(plot_df$type, levels = c("full_sample","high_ethno", "low_ethno"))
levels(plot_df$type) <- c("All Respondents\n(Prolific Sample)", "Respondents Above\nEthnocentrism Median",
                          "Respondents Below\nEthnocentrism Median")
# --- figure B2
ggplot(plot_df, aes(x = est, y = effect_type, group = country_race, color = country_race)) +
  geom_vline(xintercept = 0, size = .8, color = "black") +
  geom_linerange(aes(xmin = ci_95_low, xmax = ci_95_high), size = 1.3, color = "gray10", position = position_dodge(width = .5)) +
  geom_linerange(aes(xmin = ci_90_low, xmax = ci_90_high), size = 2.3, position = position_dodge(width = .5)) +
  geom_point(size = 5.5, color = "black", position = position_dodge(width = .5)) +
  geom_point(size = 4.5,position = position_dodge(width = .5)) +  
  facet_wrap(~type) +
  scale_color_manual(values=c("gray30", "#92978a", "#c5bec2")) +
  labs(y = NULL, x = "Estimate") +
  theme_minimal() +
  theme(legend.title = element_blank(),
        legend.text = element_text(size = 16),
        axis.text.y = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.line.y.right = element_blank(),
        axis.title.x = element_text(size = 22, margin = margin(t=15)),
        panel.spacing.x = unit(1.5, "lines"),
        panel.grid.minor.x = element_blank(),
        strip.text.x = element_text(size = 18, margin = margin (13, 0, 15, 0)),
        strip.background = element_rect(color = "black", size = 1.5, fill = "#f3f3f1"),
        legend.key = element_rect(color = NA, fill = NA),
        legend.spacing.y = unit(-8, "line"),
        legend.key.size = unit(1, "cm"))}

